2024-04-04
Cox models for time-to-event data
cph from rms to fit a Cox modelThis material is discussed in Chapters 29-31 of our Course Notes
Replicable Research and the Crisis in Science
brca <- read_csv("c22/data/brca.csv", show_col_types = FALSE) |>
mutate(across(where(is_character), as_factor),
subject = as.character(subject))
head(brca)# A tibble: 6 × 5
subject treat trial_weeks last_alive age
<chr> <fct> <dbl> <dbl> <dbl>
1 A01 S_CT 102 0 55
2 A02 S_IT 192 0 62
3 A03 S_Both 73 0 72
4 A04 S_CT 58 1 48
5 A05 S_CT 48 1 26
6 A06 S_IT 182 1 52
Data from a trial of three treatments for breast cancer
brca tibble with treat = S_CT, S_IT, S_Both and age at baselinetrial_weeks and last_alive which we used to create a survival object S.kmfit to compare the treat resultsmod_T.age) into the modelcph from the rms package to fit a Cox model that incorporates some non-linearitytrial_weeks: time in the study, in weeks, to death or censoringlast_alive: 1 if alive at last follow-up (and thus censored), 0 if deadSo last_alive = 0 if the event (death) occurs.
mod_T: Treatment alonemod_AT: Age + TreatmentCall:
coxph(formula = S ~ age + treat, data = brca)
coef exp(coef) se(coef) z p
age 0.07807 1.08119 0.03672 2.126 0.0335
treatS_IT -0.31161 0.73227 0.60936 -0.511 0.6091
treatS_Both -0.59960 0.54903 0.65741 -0.912 0.3617
Likelihood ratio test=6.99 on 3 df, p=0.07224
n= 31, number of events= 15
mod_ATtidy(mod_AT, exponentiate = TRUE, conf.int = TRUE) |>
select(term, estimate, std.error, conf.low, conf.high) |>
gt() |> fmt_number(decimals = 3) |> tab_options(table.font.size = 20)| term | estimate | std.error | conf.low | conf.high |
|---|---|---|---|---|
| age | 1.081 | 0.037 | 1.006 | 1.162 |
| treatS_IT | 0.732 | 0.609 | 0.222 | 2.417 |
| treatS_Both | 0.549 | 0.657 | 0.151 | 1.992 |
treat but Harry is one year older, the model estimates Harry will have 1.08 times the hazard of Sally (95% CI 1.01, 1.16).mod_AT| term | estimate | std.error | conf.low | conf.high |
|---|---|---|---|---|
| age | 1.081 | 0.037 | 1.006 | 1.162 |
| treatS_IT | 0.732 | 0.609 | 0.222 | 2.417 |
| treatS_Both | 0.549 | 0.657 | 0.151 | 1.992 |
S_IT and Sally receives S_CT, and they are the same age, the model estimates Cyrus will have 0.73 times the hazard of Sally (95% CI 0.22, 2.41).S_Both and Sally receives S_CT, and they are the same age, the model estimates Barry will have 0.55 times the hazard of Sally (95% CI 0.15, 1.99).n = 31, nevent = 15 for each model.
bind_rows(glance(mod_T), glance(mod_AT)) |>
mutate(model = c("mod_T", "mod_AT")) |>
select(model, p.value.log, concordance, r.squared,
max_r2 = r.squared.max, AIC, BIC) |>
gt() |> fmt_number(decimals = 3) |> tab_options(table.font.size = 20)| model | p.value.log | concordance | r.squared | max_r2 | AIC | BIC |
|---|---|---|---|---|---|---|
| mod_T | 0.416 | 0.577 | 0.055 | 0.944 | 91.773 | 93.189 |
| mod_AT | 0.072 | 0.701 | 0.202 | 0.944 | 88.536 | 90.660 |
What do the glance results indicate?
Comparing the mod_AT model with age and treatment to the mod_T model with treatment alone…
Analysis of Deviance Table
Cox model: response is S
Model 1: ~ age + treat
Model 2: ~ treat
loglik Chisq Df Pr(>|Chi|)
1 -41.268
2 -43.886 5.237 1 0.02211 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
What does this suggest? Does this match with what AIC and BIC suggested?
mod_ATcph from the rms packagerms::cph to fit a fancier AxTbrca <- read_csv("c22/data/brca.csv", show_col_types = FALSE) |>
mutate(across(where(is_character), as_factor),
subject = as.character(subject)) # reload without S
d <- datadist(brca)
options(datadist="d")
brca$S <- with(brca, Surv(trial_weeks, last_alive == 0))
cph_AxT <- cph(S ~ rcs(age, 4) + treat + age %ia% treat,
data = brca,
x = TRUE, y = TRUE, surv = TRUE)cph_AxT resultsCox Proportional Hazards Model
cph(formula = S ~ rcs(age, 4) + treat + age %ia% treat, data = brca,
x = TRUE, y = TRUE, surv = TRUE)
Model Tests Discrimination
Indexes
Obs 31 LR chi2 11.66 R2 0.332
Events 15 d.f. 7 R2(7,31) 0.140
Center 19.2233 Pr(> chi2) 0.1123 R2(7,15) 0.267
Score chi2 11.89 Dxy 0.488
Pr(> chi2) 0.1042
Coef S.E. Wald Z Pr(>|Z|)
age 0.4016 0.2610 1.54 0.1239
age' -1.2521 0.7528 -1.66 0.0963
age'' 2.7316 1.5490 1.76 0.0778
treat=S_IT 5.0537 6.1625 0.82 0.4122
treat=S_Both 4.9327 6.6650 0.74 0.4592
age * treat=S_IT -0.1011 0.1072 -0.94 0.3455
age * treat=S_Both -0.1006 0.1157 -0.87 0.3846
Effects Response : S
Factor Low High Diff. Effect S.E. Lower 0.95 Upper 0.95
age 49.5 61.5 12 0.48200 0.99998 -1.47790 2.4419
Hazard Ratio 49.5 61.5 12 1.61930 NA 0.22811 11.4950
treat - S_IT:S_CT 1.0 2.0 NA -0.40504 0.78888 -1.95120 1.1411
Hazard Ratio 1.0 2.0 NA 0.66695 NA 0.14210 3.1303
treat - S_Both:S_CT 1.0 3.0 NA -0.49745 0.80805 -2.08120 1.0863
Hazard Ratio 1.0 3.0 NA 0.60808 NA 0.12478 2.9633
Adjusted to: age=54 treat=S_CT
Divergence or singularity in 12 samples
index.orig training test optimism index.corrected n
Dxy 0.4883 0.6194 0.3722 0.2472 0.2411 28
R2 0.3320 0.5047 0.1632 0.3415 -0.0096 28
Slope 1.0000 1.0000 0.3040 0.6960 0.3040 28
D 0.1191 0.2346 0.0482 0.1864 -0.0673 28
U -0.0223 -0.0220 1.7430 -1.7650 1.7427 28
Q 0.1414 0.2566 -1.6949 1.9514 -1.8100 28
g 1.9803 9.9533 1.0352 8.9181 -6.9377 28
cph_AxT model Wald Statistics Response: S
Factor Chi-Square d.f. P
age (Factor+Higher Order Factors) 7.71 5 0.1727
All Interactions 0.96 2 0.6175
Nonlinear 3.73 2 0.1548
treat (Factor+Higher Order Factors) 2.58 4 0.6297
All Interactions 0.96 2 0.6175
age * treat (Factor+Higher Order Factors) 0.96 2 0.6175
TOTAL NONLINEAR + INTERACTION 3.74 4 0.4423
TOTAL 8.55 7 0.2868
survplot in rms for age comparisonsurvplot for treat comparisonage effect in cph_AxTtreat effect in cph_AxTcph_AxT nomogramSuppose I want to show 4-year survival rates at the bottom of the nomogram. 4 years is 208 weeks, which is the unit of time the model works with, so we have…
cph_AxT nomogramsurvminer has a function called ggcoxdiagnostics() which plots different types of residuals as a function of time, linear predictor or observation id.
type parameter, with options…type = c("martingale", "deviance", "score", "schoenfeld",
"dfbeta", "dfbetas", "scaledsch", "partial")
survminerOur department teaches an entire course on this subject every Spring (PQHS 435).
We can make acceptance of uncertainty more natural to our thinking by accompanying every point estimate in our research with a measure of its uncertainty such as a standard error or interval estimate. Reporting and interpreting point and interval estimates should be routine.
How will accepting uncertainty change anything? To begin, it will prompt us to seek better measures, more sensitive designs, and larger samples, all of which increase the rigor of research.
It also helps us be modest … [and] leads us to be thoughtful.
The nexus of openness and modesty is to report everything while at the same time not concluding anything from a single study with unwarranted certainty. Because of the strong desire to inform and be informed, there is a relentless demand to state results with certainty. Again, accept uncertainty and embrace variation in associations and effects, because they are always there, like it or not.
Understand that expressions of uncertainty are themselves uncertain. Accept that one study is rarely definitive, so encourage, sponsor, conduct, and publish replication studies.
Be modest by encouraging others to reproduce your work. Of course, for it to be reproduced readily, you will necessarily have been thoughtful in conducting the research and open in presenting it.
Switch from reliance on statistical or practical significance to the more stringent statistical criterion of practical benefit for:
Require that applied research reveal the actual unadjusted means/medians of results for all groups and subgroups, and that review panels take such data into account, as opposed to only reporting relative differences between adjusted means/medians.
Using a tidymodels approach to fit linear models.
432 Class 22 | 2024-04-04 | https://thomaselove.github.io/432-2024/